This document provides the code for the analyses reported in:

Cosme et al. (Preprint) Testing the adolescent social reorientation model using hierarchical growth curve modeling with parcellated fMRI data

This script reproduces the analyses reported in the main manuscript.

load packages

library(tidyverse)
library(knitr)
library(lme4)
library(lmerTest)

define color palettes and labels

pal_self_other = c("#FFA90A", "#247BA0")
pal_social_academic = c("#63647E", "#F25F5C")
pal_wave = c("#693668", "#A74482", "#F84AA7")
pal_label = c("#47A8BD", "#DBC057", "#FF3366")
pal_gender = c("#70c1b3","#247BA0")

parcel_labeller = labeller(label = c('social' = 'social parcels', 'other' = 'control parcels', 'self' = 'self parcels'),
                           domain = c('social' = 'social domain', 'academic' = 'academic domain'),
                           wave = c("t1" = "wave 1", "t2" = "wave 2", "t3" = "wave 3"))

label_df = expand.grid(label = c("social", "self", "other"),
              target = c("self", "other"),
              domain = c("social", "academic"),
              age = 13,
              expected_avg = 1,
              expected_diff = 1)

dcbw = theme_classic() +
  theme(text = element_text(size = 14, family = "Futura Medium", color = "black"),
        panel.background = element_blank(),
        plot.background = element_blank(),
        strip.background = element_blank(),
        strip.text = element_text(size = 14),
        legend.background = element_rect(fill = NA, color = NA),
        axis.line = element_line(color = "black"),
        axis.text = element_text(color = "black"),
        panel.grid.minor = element_blank())

load MRI data

  • exclude participants with >20% volumes with motion artifacts or low quality FX data
  • standardize betas within parcel (across participants and time)
  • recode standardized betas +/- 3 SDs from the mean as NA (0.8% of all observations)
# define parcels
self_parcels = c(5, 17, 28, 47, 52, 66, 116, 147, 152, 156, 184, 198, 207, 225, 249, 292, 309, 354, 380)
social_parcels = c(18, 23, 29, 49, 54, 59, 62, 63, 67, 76, 111, 114, 139, 143, 146, 150, 178, 179, 189, 203, 212, 224, 229, 236, 238, 239, 245, 250, 259, 266, 271, 301, 305, 310, 322, 324, 328, 331, 333, 342, 343, 350, 374, 391)

# mri exclusions
mri_exclusions = c('s002_t1', 's004_t1', 's008_t1', 's011_t1', 's017_t1', 
                   's026_t1', 's033_t2', 's034_t1', 's041_t1', 's044_t1', 
                   's047_t1', 's051_t1', 's054_t1', 's057_t1', 's059_t1', 
                   's061_t1', 's063_t1', 's070_t2', 's074_t1', 's074_t2', 
                   's078_t1', 's084_t1', 's090_t2', 's090_t3', 's094_t1', 
                   's094_t2', 's096_t1') 

# load and tidy parcel data
parcellations = read_csv('../data/fxParcellations.csv') %>%
  mutate(label = ifelse(parcellation %in% self_parcels, 'self',
                 ifelse(parcellation %in% social_parcels, 'social', 'other'))) %>%
  mutate(wave = paste0("t", as.numeric(c(`10` = 1, `13` = 2, `16` = 3)[as.character(age)]))) %>%
  select(-age) %>%
  unite(sub_wave, c(subjectID, wave), remove = FALSE) %>%
  group_by(parcellation) %>%
  mutate(inclusion = ifelse(sub_wave %in% mri_exclusions, "excluded from MRI", "completed MRI"),
         beta = ifelse(sub_wave %in% mri_exclusions, NA, beta),
         sd = ifelse(sub_wave %in% mri_exclusions, NA, sd),
         beta_std = scale(beta, center = FALSE, scale = TRUE),
         mean_beta_std = mean(beta_std, na.rm = TRUE)) %>%
    select(-sub_wave) %>%
  ungroup()

# exclude parameter estimates 3 SD from the mean
parcellations_ex = parcellations %>%
  mutate(beta_std = ifelse(beta_std > mean_beta_std + 3 | beta_std < mean_beta_std - 3, NA, beta_std))

load demographic data and merge

# demographics
demo = read.csv("../data/SFIC_age.pds.gender.csv") %>%
  rename("subjectID" = SID,
         "wave" = wavenum,
         "gender" = Gender) %>%
  mutate(subjectID = sprintf("s%03d", subjectID),
         wave = paste0("t", wave),
         age_c = age - 13,
         age_c2 = age_c^2,
         pdss_c = pdss - 3,
         pdss_c2 = pdss_c^2)

# merge data
merged = parcellations_ex %>%
  full_join(., demo, by = c("subjectID", "wave")) %>%
  mutate(inclusion = ifelse(is.na(inclusion), "didn't complete MRI", inclusion)) %>%
  filter(!(subjectID == "s086" & wave == "t3")) #no MRI, task, or self-report data was collected

# subset data for modeling
neuro_model_data = merged %>%
  filter(!is.na(beta)) %>%
  select(subjectID, wave, age, age_c, age_c2, target, domain, parcellation, label, beta, beta_std)

neuro_model_data_ex = neuro_model_data  %>%
  na.omit()

# dummy code target and domain
neuro_model_data_dummy = neuro_model_data_ex %>%
  mutate(target = ifelse(target == "self", .5, -.5),
         domain = ifelse(domain == "social", .5, -.5))

sample descriptives

full sample

merged %>% 
  select(subjectID, wave, age, gender) %>%
  unique() %>%
  group_by(wave) %>%
  summarize(n = n(),
            meanAge = mean(age, na.rm = TRUE),
            sdAge = sd(age, na.rm = TRUE)) %>%
  kable(format = 'pandoc', digits = 2)
wave n meanAge sdAge
t1 90 10.07 0.31
t2 57 13.04 0.31
t3 44 16.46 0.58
merged %>% 
  select(subjectID, wave, age, gender) %>%
  unique() %>%
  group_by(wave, gender) %>%
  summarize(n = n()) %>%
  spread(gender, n) %>%
  kable(format = 'pandoc', digits = 2)
wave F M
t1 45 45
t2 30 27
t3 25 19

broken down by inclusion

age_table = merged %>% 
  select(subjectID, inclusion, wave, age, pdss, gender) %>%
  unique() %>%
  group_by(inclusion, wave) %>%
  summarize(N = n(),
            Mage = mean(age, na.rm = TRUE),
            SDage = sd(age, na.rm = TRUE),
            Mpdss = mean(pdss, na.rm = TRUE),
            SDpdss = sd(pdss, na.rm = TRUE))

gender_table = merged %>% 
  select(subjectID, inclusion, wave, age, pds, gender) %>%
  unique() %>%
  group_by(inclusion, wave, gender) %>%
  summarize(N = n()) %>%
  spread(gender, N) %>%
  rename("Female" = F,
         "Male" = M)

merged %>% 
  filter(!inclusion == "didn't complete MRI") %>%
  mutate(inclusion = ifelse(grepl("excluded", inclusion), "excluded", "included")) %>%
  select(subjectID, wave, gender, age, inclusion) %>%
  unique() %>%
  filter(inclusion == "included") %>%
  group_by(subjectID) %>%
  summarize(`number of waves included` = n()) %>%
  group_by(`number of waves included`) %>%
  summarize(n = n()) %>%
  kable(format = 'pandoc')
number of waves included n
1 35
2 17
3 22
age_table %>%
  left_join(gender_table) %>%
  extract(wave, "Wave", "t([1-3]{1})") %>%
  arrange(Wave) %>%
  select(Wave, inclusion, N, Mage, SDage, Mpdss, SDpdss, Female, Male) %>%
  rename("MRI status" = inclusion) %>%
  kable(format = 'pandoc', digits = 2)
Wave MRI status N Mage SDage Mpdss SDpdss Female Male
1 completed MRI 57 10.08 0.32 2.16 0.74 33 24
1 didn’t complete MRI 12 10.00 0.25 2.00 0.77 4 8
1 excluded from MRI 21 10.07 0.33 1.83 0.58 8 13
2 completed MRI 44 13.06 0.33 3.57 1.05 25 19
2 didn’t complete MRI 8 13.07 0.20 3.69 0.65 4 4
2 excluded from MRI 5 12.79 0.28 2.62 0.85 1 4
3 completed MRI 34 16.33 0.46 4.75 0.39 19 15
3 didn’t complete MRI 9 17.06 0.60 4.78 0.67 6 3
3 excluded from MRI 1 15.72 NA 3.50 NA NA 1

plot

full sample

merged %>% 
  select(subjectID, wave, gender, age, inclusion) %>%
  group_by(subjectID) %>%
  unique() %>%
  mutate(n.waves = n()) %>%
  ungroup() %>%
  arrange(n.waves, gender, age, inclusion) %>%
  mutate(order = row_number()) %>%
  ggplot(aes(age, reorder(subjectID, -order), color = gender)) +
    geom_line(size = .3) + 
    geom_point(aes(shape = inclusion), size = 2, fill = "white") +
    scale_color_manual(values = pal_gender, name = "", labels = c("female", "male")) +
    scale_shape_manual(name = "", values = c(19, 17, 21)) +
    labs(x = "\nage", y = "participant\n") +
    dcbw +
    theme(axis.line = element_line(colour = "black"),
          axis.text.y = element_blank(),
          axis.ticks.y = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          legend.position = c(.75,.92),
          legend.spacing.y = unit(.01, 'cm'),
          legend.margin = unit(0, "cm"))

MRI only

merged %>% 
  filter(!inclusion == "didn't complete MRI") %>%
  mutate(inclusion = ifelse(grepl("excluded", inclusion), "excluded", "included")) %>%
  select(subjectID, wave, gender, age, inclusion) %>%
  group_by(subjectID) %>%
  unique() %>%
  mutate(n.waves = n()) %>%
  ungroup() %>%
  arrange(n.waves, gender, age, inclusion) %>%
  mutate(order = row_number()) %>%
  ggplot(aes(age, reorder(subjectID, -order), color = gender)) +
    geom_line(size = .3) + 
    geom_point(aes(shape = inclusion), size = 2, fill = "white") +
    scale_color_manual(values = pal_gender, name = "", labels = c("female", "male")) +
    scale_shape_manual(name = "", values = c(21, 19)) +
    labs(x = "\nage", y = "participant\n") +
    dcbw +
    theme(axis.line = element_line(colour = "black"),
          axis.text.y = element_blank(),
          axis.ticks.y = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          legend.position = c(.75,.92),
          legend.spacing.y = unit(.01, 'cm'),
          legend.margin = unit(0, "cm"))

mri descriptives

# collapsed across wave
neuro_model_data_dummy %>%
  mutate(target = ifelse(target == .5, "self", "other"),
         domain = ifelse(domain == .5, "social", "academic")) %>%
  group_by(domain, label) %>%
  summarize(M = mean(beta_std, na.rm = TRUE),
            SD = sd(beta_std, na.rm = TRUE)) %>%
  gather(var, value, M, SD) %>%
  unite(label, label, var, sep = " ") %>%
  spread(label, value) %>%
  kable(format = "pandoc", digits = 2)
domain other M other SD self M self SD social M social SD
academic 0.01 0.92 -0.08 0.91 0.14 0.92
social -0.01 0.95 0.02 0.92 0.27 0.93
neuro_model_data_dummy %>%
  mutate(target = ifelse(target == .5, "self", "other"),
         domain = ifelse(domain == .5, "social", "academic")) %>%
  group_by(target, label) %>%
  summarize(M = mean(beta_std, na.rm = TRUE),
            SD = sd(beta_std, na.rm = TRUE)) %>%
  gather(var, value, M, SD) %>%
  unite(label, label, var, sep = " ") %>%
  spread(label, value) %>%
  kable(format = "pandoc", digits = 2)
target other M other SD self M self SD social M social SD
other 0.00 0.94 -0.11 0.91 0.22 0.93
self -0.01 0.94 0.05 0.92 0.19 0.93
# as a function of wave
neuro_model_data_dummy %>%
  mutate(target = ifelse(target == .5, "self", "other"),
         domain = ifelse(domain == .5, "social", "academic")) %>%
  group_by(domain, label, wave) %>%
  summarize(M = mean(beta_std, na.rm = TRUE),
            SD = sd(beta_std, na.rm = TRUE)) %>%
  gather(var, value, M, SD) %>%
  unite(label, label, var, sep = " ") %>%
  spread(label, value) %>%
  kable(format = "pandoc", digits = 2)
domain wave other M other SD self M self SD social M social SD
academic t1 0.00 0.92 -0.13 0.95 0.06 0.91
academic t2 0.00 0.96 -0.04 0.92 0.19 0.95
academic t3 0.03 0.89 -0.04 0.82 0.21 0.88
social t1 -0.02 0.94 0.00 0.91 0.22 0.91
social t2 0.00 0.96 0.00 0.95 0.30 0.98
social t3 -0.01 0.96 0.05 0.89 0.33 0.90
neuro_model_data_dummy %>%
  mutate(target = ifelse(target == .5, "self", "other"),
         domain = ifelse(domain == .5, "social", "academic")) %>%
  group_by(target, label, wave) %>%
  summarize(M = mean(beta_std, na.rm = TRUE),
            SD = sd(beta_std, na.rm = TRUE)) %>%
  gather(var, value, M, SD) %>%
  unite(label, label, var, sep = " ") %>%
  spread(label, value) %>%
  kable(format = "pandoc", digits = 2)
target wave other M other SD self M self SD social M social SD
other t1 0.00 0.92 -0.10 0.92 0.16 0.90
other t2 0.00 0.97 -0.14 0.94 0.27 0.98
other t3 0.02 0.92 -0.10 0.85 0.27 0.88
self t1 -0.02 0.94 -0.03 0.94 0.11 0.92
self t2 0.00 0.94 0.11 0.92 0.22 0.95
self t3 0.01 0.93 0.11 0.86 0.27 0.90

visualize raw

Visualize the developmental trajectories using the raw data

hypothesis 0

domain_parc_plot_raw = neuro_model_data %>%
  distinct(parcellation, target, domain, age, label, .keep_all = T) %>%
  group_by(subjectID, wave, label, domain, parcellation) %>%
  mutate(beta_std_avg = mean(beta_std, na.rm = TRUE)) %>%
  ggplot(aes(x = age, y = beta_std_avg, color = domain)) + 
  geom_smooth(aes(group = interaction(parcellation, domain), size = label),
              method = "lm", formula = y ~ poly(x, 2), se = FALSE, show.legend = FALSE) +
  geom_smooth(aes(fill = domain), method = "lm", formula = y ~ poly(x, 2), size = 1.5) +
  scale_color_manual(values = pal_social_academic) +
  scale_fill_manual(values = pal_social_academic) +
  scale_size_manual(values = c(.05, .1, .1)) +
  scale_x_continuous(breaks = c(10, 13, 16)) + 
  scale_y_continuous(breaks = c(-1, 0, 1)) + 
  coord_cartesian(ylim = c(-1.2, 1.2)) +
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = '\nage', y = 'mean standardized parameter estimate\n', color = '', fill = '') +
  theme_minimal(base_size = 12) +
  dcbw +
  theme(legend.position = c(.85, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(0, "cm"))

target_parc_plot_raw = neuro_model_data %>%
  distinct(parcellation, target, domain, age, label, .keep_all = T) %>%
  group_by(subjectID, wave, label, target, parcellation) %>%
  mutate(beta_std_avg = mean(beta_std, na.rm = TRUE)) %>%
  ggplot(aes(x = age, y = beta_std_avg, color = target)) + 
  geom_smooth(aes(group = interaction(parcellation, target), size = label),
              method = "lm", formula = y ~ poly(x, 2), se = FALSE, show.legend = FALSE) +
  geom_smooth(aes(fill = target), method = "lm", formula = y ~ poly(x, 2), size = 1.5) +
  scale_color_manual(values = pal_self_other) +
  scale_fill_manual(values = pal_self_other) +
  scale_size_manual(values = c(.05, .1, .1)) +
  scale_x_continuous(breaks = c(10, 13, 16)) + 
  scale_y_continuous(breaks = c(-1, 0, 1)) + 
  coord_cartesian(ylim = c(-1.2, 1.2)) +
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = '\nage', y = 'mean standardized parameter estimate\n', color = '', fill = '') +
  theme_minimal(base_size = 12) +
  dcbw +
  theme(legend.position = c(.85, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(0, "cm"))

(h0_raw = cowplot::plot_grid(domain_parc_plot_raw, target_parc_plot_raw,
                                labels = c('A', 'B'), ncol = 2,
                                rel_widths = c(1, 1)))

hypothesis 1

domain_plot_raw = neuro_model_data %>%
  distinct(parcellation, target, domain, age, label, .keep_all = T) %>%
  group_by(subjectID, wave, label, domain, parcellation) %>%
  mutate(beta_std_avg = mean(beta_std, na.rm = TRUE)) %>%
  ggplot(aes(x = age, y = beta_std_avg, group = domain, color = domain)) +
  geom_smooth(aes(fill = domain), method = 'lm', formula = y ~ poly(x,2), alpha = .2) + 
  scale_color_manual(values = pal_social_academic) +
  scale_fill_manual(values = pal_social_academic) +
  scale_x_continuous(breaks = c(10, 13, 16)) + 
  scale_y_continuous(breaks = round(seq(-.4, .6, .2), 1)) +
  coord_cartesian(ylim = c(-.45, .65)) +
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = '\nage', y = 'mean standardized parameter estimate\n', color = '', fill = '') +
  dcbw +
  theme(legend.position = "none")

soc_acad_plot_raw = neuro_model_data %>%
  group_by(subjectID, age, label, domain, parcellation) %>%
  mutate(beta_std_avg = mean(beta_std, na.rm = TRUE)) %>%
  select(subjectID, age, label, domain, beta_std_avg) %>%
  unique() %>%
  spread(domain, beta_std_avg) %>%
  mutate(avg_diff = social - academic) %>%
  distinct(parcellation, age, label, .keep_all = T) %>%
  ggplot(aes(x = age, y = avg_diff)) +
  geom_smooth(aes(group = parcellation, size = label), method = "lm", formula = y ~ poly(x, 2),
              se = FALSE, color = "grey50") +
  geom_smooth(method = 'lm', formula = y ~ poly(x,2), 
              se = TRUE, color = pal_social_academic[2], fill = pal_social_academic[2], size = 1.5) + 
  scale_size_manual(values = c(.03, .1, .1)) +
  scale_y_continuous(breaks = round(seq(-.4, .6, .2), 1)) +
  coord_cartesian(ylim = c(-.45, .65)) +
  scale_x_continuous(breaks = c(10, 13, 16)) + 
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\nage", y = "mean standardized BOLD signal value\n",  color = '', fill = '') +
  dcbw +
  theme(legend.position = "none")

(h1_raw = cowplot::plot_grid(domain_plot_raw, soc_acad_plot_raw,
                                labels = c('A', 'B'), ncol = 2,
                                rel_widths = c(1, 1)))

hypothesis 2

target_plot_raw = neuro_model_data %>%
  distinct(parcellation, target, domain, age, label, .keep_all = T) %>%
  group_by(subjectID, wave, label, target, parcellation) %>%
  mutate(beta_std_avg = mean(beta_std, na.rm = TRUE)) %>%
  ggplot(aes(x = age, y = beta_std_avg, group = target, color = target)) +
  geom_smooth(aes(fill = target), method = 'lm', formula = y ~ poly(x,2), alpha = .2) + 
  scale_color_manual(values = pal_self_other) +
  scale_fill_manual(values = pal_self_other) +
  scale_x_continuous(breaks = c(10, 13, 16)) + 
  scale_y_continuous(breaks = round(seq(-.4, .4, .2), 1)) +
  coord_cartesian(ylim = c(-.4, .5)) +
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = '\nage', y = 'mean standardized parameter estimate\n', color = '', fill = '') +
  dcbw +
  theme(legend.position = "none")

self_other_plot_raw = neuro_model_data %>%
  group_by(subjectID, age, label, target, parcellation) %>%
  mutate(beta_std_avf = mean(beta_std, na.rm = TRUE)) %>%
  select(subjectID, age, label, target, beta_std_avf) %>%
  unique() %>%
  spread(target, beta_std_avf) %>%
  mutate(avg_diff = self - other) %>%
  distinct(parcellation, age, label, .keep_all = T) %>%
  ggplot(aes(x = age, y = avg_diff)) +
  geom_smooth(aes(group = parcellation, size = label), method = "lm", formula = y ~ poly(x, 2),
              se = FALSE, color = "grey50") +
  geom_smooth(method = 'lm', formula = y ~ poly(x,2), 
              se = TRUE, color = pal_self_other[2], fill = pal_self_other[2], size = 1.5) + 
  scale_size_manual(values = c(.03, .1, .1)) +
  scale_y_continuous(breaks = round(seq(-.4, .4, .2), 1)) +
  coord_cartesian(ylim = c(-.4, .5)) +
  scale_x_continuous(breaks = c(10, 13, 16)) + 
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\nage", y = "mean standardized BOLD signal value\n",  color = '', fill = '') +
  dcbw +
  theme(legend.position = "none")

(h2_raw = cowplot::plot_grid(target_plot_raw, self_other_plot_raw,
                                labels = c('A', 'B'), ncol = 2,
                                rel_widths = c(1, 1)))

hypothesis 3

int_plot_raw = neuro_model_data %>%
  distinct(parcellation, target, domain, age, beta_std, label, .keep_all = T) %>%
  ggplot(aes(x = age, y = beta_std, group = interaction(target, domain), color = domain, linetype = target)) +
  geom_smooth(aes(fill = domain), method = 'lm', formula = y ~ poly(x,2), alpha = .2) + 
  scale_y_continuous(breaks = c(-.4, -.2, 0, .2, .4)) +
  coord_cartesian(ylim = c(-.4, .45)) +
  scale_color_manual(values = pal_social_academic) +
  scale_fill_manual(values = pal_social_academic) +
  scale_linetype_manual(name = "", values = c("dotted", "solid")) +
  scale_x_continuous(breaks = c(10, 13, 16)) + 
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = '\nage', y = 'mean standardized parameter estimate\n', color = '', linetype = '', fill = '') +
  guides(linetype = guide_legend(override.aes = list(color = "black"))) +
  dcbw +
  theme(legend.position = c(.75, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(.5, "cm"),
        legend.direction = "vertical",
        legend.box = "horizontal")

soc_acad_int_plot_raw = neuro_model_data %>%
  group_by(subjectID, age, label, target, domain, parcellation) %>%
  mutate(beta_std_avg = mean(beta_std, na.rm = TRUE)) %>%
  select(subjectID, age, label, target, domain, beta_std_avg) %>%
  unique() %>%
  spread(domain, beta_std_avg) %>%
  mutate(avg_diff = social - academic) %>%
  distinct(parcellation, age, label, .keep_all = T) %>%
  ggplot(aes(x = age, y = avg_diff, color = target)) +
  geom_smooth(aes(group = interaction(parcellation, target), size = label),
              method = "lm", formula = y ~ poly(x, 2), se = FALSE, show.legend = FALSE) +
  geom_smooth(aes(fill = target), method = 'lm', formula = y ~ poly(x,2), 
              se = TRUE, size = 1.5) + 
  scale_color_manual(values = pal_self_other) +
  scale_fill_manual(values = pal_self_other) +
  scale_size_manual(values = c(.03, .1, .1)) +
  scale_y_continuous(breaks = round(seq(-.6, .6, .2), 1)) +
  coord_cartesian(ylim = c(-.65, .65)) +
  scale_x_continuous(breaks = c(10, 13, 16)) + 
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\nage", y = "mean standardized BOLD signal value\n",  color = '', fill = '') +
  dcbw +
  theme(legend.position = c(.85, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(0, "cm"))

self_other_int_plot_raw = neuro_model_data %>%
  group_by(subjectID, age, label, domain, target, parcellation) %>%
  mutate(beta_std_avg = mean(beta_std, na.rm = TRUE)) %>%
  select(subjectID, age, label, domain, target, beta_std_avg) %>%
  unique() %>%
  spread(target, beta_std_avg) %>%
  mutate(avg_diff = self - other) %>%
  distinct(parcellation, age, label, .keep_all = T) %>%
  ggplot(aes(x = age, y = avg_diff, color = domain)) +
  geom_smooth(aes(group = interaction(parcellation, domain), size = label),
              method = "lm", formula = y ~ poly(x, 2), se = FALSE, show.legend = FALSE) +
  geom_smooth(aes(fill = domain), method = 'lm', formula = y ~ poly(x,2), 
              se = TRUE, size = 1.5) + 
  scale_color_manual(values = pal_social_academic) +
  scale_fill_manual(values = pal_social_academic) +
  scale_size_manual(values = c(.03, .1, .1)) +
  scale_y_continuous(breaks = round(seq(-.6, .6, .2), 1)) +
  coord_cartesian(ylim = c(-.65, .65)) +
  scale_x_continuous(breaks = c(10, 13, 16)) + 
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\nage", y = "mean standardized BOLD signal value\n",  color = '', fill = '') +
  dcbw +
  theme(legend.position = c(.85, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(0, "cm"))

(h3_raw = cowplot::plot_grid(int_plot_raw, soc_acad_int_plot_raw, self_other_int_plot_raw,
                       labels = c('A', 'B', 'C'), ncol = 3))

run models

Estimate full models

domain x age

# specify model
model_domain_equation = formula(beta_std ~ 1 + age_c*domain*label +
                                  age_c2*domain*label +
                                  (1 + domain + age_c | subjectID) +
                                  (1 + domain + age_c + age_c2 | parcellation))

model_domain_formula = lFormula(model_domain_equation, data = neuro_model_data_dummy)

# calculate max number of iterations
model_domain_numFx = length(dimnames(model_domain_formula$X)[[2]])

model_domain_numRx = sum(as.numeric(lapply(model_domain_formula$reTrms$cnms, function(x) {
  l <- length(x)
  (l*(l - 1)) / 2 + l
})))

model_domain_maxfun = 50*(model_domain_numFx + model_domain_numRx + 1)^2

# run or load the model
if (file.exists("../data/model_domain.RDS")) {
  model_domain = readRDS("../data/model_domain.RDS")
} else {
  model_domain = lmer(model_domain_equation, data = neuro_model_data_dummy, REML = F, #Use ML since we want to compare random effects
                      verbose = 2,
                      control = lmerControl(optCtrl = list(maxfun = model_domain_maxfun), optimizer = "bobyqa", calc.derivs = FALSE))
  saveRDS(model_domain, "../data/model_domain.RDS")
}

domain x target x age

# specify model
model_target_equation = formula(beta_std ~ 1 + age_c*target*domain*label +
                                  age_c2*target*domain*label +
                                  (1 + age_c*target*domain | subjectID) +
                                  (1 + age_c*target*domain + age_c2*target*domain | parcellation))

# calculate max number of iterations
model_target_formula = lFormula(model_target_equation, data = neuro_model_data_dummy)

model_target_numFx = length(dimnames(model_target_formula$X)[[2]])

model_target_numRx = sum(as.numeric(lapply(model_target_formula$reTrms$cnms, function(x) {
  l <- length(x)
  (l*(l - 1)) / 2 + l
})))

model_target_maxfun = 10*(model_target_numFx + model_target_numRx + 1)^2

# run or load the model
if (file.exists("../data/model_target.RDS")) {
  model_target = readRDS("../data/model_target.RDS")
} else {
  model_target = lmer(model_target_equation, data = neuro_model_data_dummy, REML = F, #Use ML since we want to compare random effects
                      verbose = 2,
                      control = lmerControl(optCtrl = list(maxfun = model_target_maxfun), optimizer = "bobyqa", calc.derivs = FALSE))
  saveRDS(model_target, "../data/model_target.RDS")
}

compare models

anova(model_domain, model_target) %>%
  as.data.frame() %>%
  select(npar, AIC, BIC, Chisq, Df, `Pr(>Chisq)`) %>%
  rownames_to_column() %>%
  rename("model" = rowname,
         "model df" = Df,
         "X2" = Chisq,
         "X2 df" = Df,
         "p" = `Pr(>Chisq)`) %>%
  mutate(AIC = round(AIC, 2),
         BIC = round(BIC, 2),
         X2 = round(X2, 2),
         p = round(p, 3),
         p = ifelse(p == 0, "< .001", gsub("0.(.*)", ".\\1", p))) %>%
  mutate_if(is.numeric, funs(ifelse(is.na(.), "--", .))) %>%
  mutate_if(is.character, funs(ifelse(is.na(.), "--", .))) %>%
  kable(format = "pandoc")
model npar AIC BIC X2 X2 df p
model_domain 35 462224.6 462579.8 – – –
model_target 151 461839.9 463372.1 616.74 116 < .001

summarize best fitting model

model_target %>%
  broom.mixed::tidy(effects = c("ran_pars", "fixed"), conf.int = TRUE) %>%
  filter(effect == "fixed") %>%
  select(-group) %>%
  rename("b" = estimate,
         "SE" = std.error,
         "t" = statistic,
         "p" = p.value) %>%
  mutate(p = round(p, 3),
         p = ifelse(p == 0, "< .001", gsub("0.(.*)", ".\\1", sprintf("%.3f", p))),
         term = gsub("\\(Intercept\\)", "Intercept (age 13, label (control))", term),
         term = gsub("target", "Target", term),
         term = gsub("domain", "Domain", term),
         term = gsub("labelself", "Label (self)", term),
         term = gsub("labelsocial", "Label (social)", term),
         term = gsub("age_c", "Age", term),
         term = gsub(":", " x ", term),
         term = gsub("sd__", "", term),
         term = gsub("Observation", "observation", term),
         effect = gsub("ran_pars", "random", effect),
         `b [95% CI]` = ifelse(effect == "fixed", sprintf("%.3f [%.3f, %.3f]", b, conf.low, conf.high), "--")) %>%
  mutate_if(is.numeric, round, 3) %>%
  mutate_if(is.numeric, funs(ifelse(is.na(.), "--", .))) %>%
  mutate_if(is.character, funs(ifelse(is.na(.), "--", .))) %>%
  select(term, `b [95% CI]`, SE, t, df, p) %>%
  kable(format = "pandoc")
term b [95% CI] SE t df p
Intercept (age 13, label (control)) -0.004 [-0.058, 0.050] 0.028 -0.142 380.663 .887
Age 0.003 [-0.002, 0.009] 0.003 1.169 146.119 .244
Target 0.004 [-0.015, 0.023] 0.010 0.399 188.970 .690
Domain 0.013 [-0.012, 0.039] 0.013 1.036 139.124 .302
Label (self) -0.033 [-0.247, 0.180] 0.109 -0.307 352.122 .759
Label (social) 0.242 [0.096, 0.388] 0.074 3.253 352.073 .001
Age2 0.000 [-0.001, 0.001] 0.001 0.217 577.245 .828
Age x Target 0.002 [-0.002, 0.007] 0.002 0.995 75.175 .323
Age x Domain 0.000 [-0.007, 0.008] 0.004 0.133 58.010 .895
Target x Domain -0.006 [-0.043, 0.030] 0.019 -0.325 137.302 .746
Age x Label (self) 0.007 [-0.008, 0.022] 0.008 0.914 352.492 .361
Age x Label (social) 0.018 [0.008, 0.029] 0.005 3.470 352.023 .001
Target x Label (self) 0.255 [0.192, 0.318] 0.032 7.920 475.077 < .001
Target x Label (social) -0.051 [-0.094, -0.008] 0.022 -2.319 472.033 .021
Domain x Label (self) 0.066 [-0.006, 0.138] 0.037 1.804 447.869 .072
Domain x Label (social) 0.118 [0.069, 0.167] 0.025 4.713 445.664 < .001
Target x Age2 -0.002 [-0.004, 0.000] 0.001 -1.813 959.625 .070
Domain x Age2 -0.003 [-0.005, -0.001] 0.001 -3.004 1434.413 .003
Label (self) x Age2 0.001 [-0.003, 0.006] 0.002 0.577 366.801 .564
Label (social) x Age2 -0.004 [-0.007, -0.001] 0.002 -2.821 364.916 .005
Age x Target x Domain 0.003 [-0.008, 0.013] 0.005 0.525 54.843 .602
Age x Target x Label (self) 0.026 [0.013, 0.039] 0.007 3.843 3904.827 < .001
Age x Target x Label (social) 0.006 [-0.003, 0.015] 0.005 1.365 3877.736 .172
Age x Domain x Label (self) -0.006 [-0.020, 0.008] 0.007 -0.861 1437.557 .389
Age x Domain x Label (social) -0.005 [-0.014, 0.005] 0.005 -0.955 1428.293 .340
Target x Domain x Label (self) -0.022 [-0.134, 0.090] 0.057 -0.385 3796.005 .700
Target x Domain x Label (social) 0.098 [0.021, 0.174] 0.039 2.499 3765.899 .013
Target x Domain x Age2 -0.000 [-0.004, 0.004] 0.002 -0.067 1743.880 .946
Target x Label (self) x Age2 -0.011 [-0.018, -0.004] 0.004 -3.014 1444.407 .003
Target x Label (social) x Age2 0.004 [-0.000, 0.009] 0.002 1.764 1432.563 .078
Domain x Label (self) x Age2 0.006 [-0.000, 0.013] 0.003 1.876 3673.416 .061
Domain x Label (social) x Age2 0.005 [0.001, 0.010] 0.002 2.327 3641.719 .020
Age x Target x Domain x Label (self) -0.006 [-0.032, 0.021] 0.013 -0.417 17643.185 .677
Age x Target x Domain x Label (social) -0.014 [-0.031, 0.004] 0.009 -1.497 17526.368 .134
Target x Domain x Label (self) x Age2 0.009 [-0.005, 0.022] 0.007 1.247 10786.627 .212
Target x Domain x Label (social) x Age2 -0.001 [-0.010, 0.008] 0.005 -0.290 10695.055 .772
# summarize random effects
print(VarCorr(model_target), comp = c("Variance","Std.Dev."), digits = 3)
##  Groups       Name                 Variance   Std.Dev. Corr                   
##  parcellation (Intercept)          0.20831548 0.45642                         
##               age_c                0.00088442 0.02974   0.46                  
##               target               0.00450060 0.06709  -0.11 -0.35            
##               domain               0.00999328 0.09997   0.23  0.30 -0.32      
##               age_c2               0.00003842 0.00620  -0.37 -0.37  0.53 -0.30
##               age_c:target         0.00002869 0.00536  -0.30  0.06  0.78 -0.47
##               age_c:domain         0.00009994 0.01000   0.75  0.58 -0.02 -0.33
##               target:domain        0.00238388 0.04882   0.31  0.71 -0.60  0.84
##               target:age_c2        0.00002116 0.00460   0.43  0.43 -0.86  0.19
##               domain:age_c2        0.00000705 0.00265   0.22 -0.28  0.00 -0.65
##               age_c:target:domain  0.00002628 0.00513   0.65 -0.10  0.38 -0.45
##               target:domain:age_c2 0.00001158 0.00340  -0.85 -0.59  0.56 -0.26
##  subjectID    (Intercept)          0.00198658 0.04457                         
##               age_c                0.00019403 0.01393   0.18                  
##               target               0.00141727 0.03765  -0.16  0.13            
##               domain               0.00473264 0.06879   0.10  0.31  0.04      
##               age_c:target         0.00015482 0.01244  -0.09 -0.07 -0.11 -0.39
##               age_c:domain         0.00054099 0.02326  -0.50  0.04  0.31  0.18
##               target:domain        0.00739223 0.08598   0.18  0.45  0.11 -0.32
##               age_c:target:domain  0.00089498 0.02992   0.06 -0.08 -0.43 -0.56
##  Residual                          0.66422200 0.81500                         
##                                           
##                                           
##                                           
##                                           
##                                           
##                                           
##   0.48                                    
##  -0.37  0.12                              
##  -0.52 -0.43  0.02                        
##  -0.82 -0.72  0.42  0.52                  
##  -0.41 -0.11  0.51 -0.52  0.38            
##  -0.17  0.10  0.72 -0.46  0.08  0.74      
##   0.69  0.56 -0.71 -0.52 -0.83 -0.31 -0.41
##                                           
##                                           
##                                           
##                                           
##                                           
##   0.11                                    
##   0.33 -0.08                              
##   0.44 -0.74  0.30                        
## 

simple slopes

Estimate simple slopes to test interactions at specific levels

run models

# self social > academic
self_social = emmeans::emtrends(model_target, pairwise ~ domain,
                  var = "age_c", at = list(target =.5, label="social"),
                  lmerTest.limit = 188577)$contrasts %>%
  data.frame() %>%
  mutate(contrast = "self social > academic",
         parcel = "social",
         age_effect = "linear") %>%
  bind_rows(emmeans::emtrends(model_target, pairwise ~ domain,
                  var = "age_c2", at = list(target =.5, label="social"),
                  lmerTest.limit = 188577)$contrasts %>%
  data.frame() %>%
  mutate(contrast = "self social > academic",
         parcel = "social",
         age_effect = "quadratic"))

# social self > other
social_self = emmeans::emtrends(model_target, pairwise ~ target,
                  var = "age_c", at = list(domain =.5, label="self"),
                  lmerTest.limit = 188577)$contrasts %>%
  data.frame() %>%
  mutate(contrast = "social self > other",
         parcel = "self",
         age_effect = "linear") %>%
  bind_rows(emmeans::emtrends(model_target, pairwise ~ target,
                  var = "age_c2", at = list(domain =.5, label="self"),
                  lmerTest.limit = 188577)$contrasts %>%
  data.frame() %>%
  mutate(contrast = "social self > other",
         parcel = "self",
         age_effect = "quadratic"))

make table

social_self %>%
  bind_rows(self_social) %>%
  select(contrast, parcel, age_effect, estimate, SE, df, t.ratio, p.value) %>%
  rename("b" = estimate,
         "t" = t.ratio,
         "p" = p.value) %>%
  mutate(b = round(b, 3) * -1, #flip signs for it's .5 - (-.5)
         SE = round(SE, 3),
         df = round(df, 2),
         t = abs(round(t, 2)),
         p = round(p, 3)) %>%
  kable(format = "pandoc")
contrast parcel age_effect b SE df t p
social self > other self linear 0.027 0.010 1664.99 2.77 0.006
social self > other self quadratic -0.008 0.005 4879.37 1.75 0.081
self social > academic social linear -0.010 0.007 910.55 1.43 0.153
self social > academic social quadratic 0.002 0.003 7651.76 0.53 0.598

visualize fitted models

Visualize the developmental trajectory using the fitted values from the domain x target x age model

reForm = as.formula("~(1 + age_c*target*domain + age_c2*target*domain | parcellation)")

neuro_plot_data = with(neuro_model_data_dummy, 
                    expand.grid(target = unique(target), 
                                domain = unique(domain),
                                parcellation = unique(parcellation),
                                age = unique(age),
                                stringsAsFactors = F)) %>%
  mutate(label = ifelse(parcellation %in% self_parcels, 'self',
                 ifelse(parcellation %in% social_parcels, 'social', 'other')),
         age_c = age - 13, 
         age_c2 = age_c^2, 
         subjectID = NA)
neuro_plot_data$expected = predict(model_target, newdata = neuro_plot_data, re.form = reForm)
neuro_plot_data$expected_mean = predict(model_target, newdata = neuro_plot_data, re.form = NA)

neuro_plot_data = neuro_plot_data %>%
  mutate(target = factor(target, levels = c(-.5, .5), labels = c("other", "self")),
         domain = factor(domain, levels = c(-.5, .5), labels = c("academic", "social")))

hypothesis 0

domain_parc_plot = neuro_plot_data %>%
  distinct(parcellation, target, domain, age, label, .keep_all = T) %>%
  group_by(subjectID, age, label, domain, parcellation) %>%
  mutate(expected_avg = mean(expected, na.rm = TRUE)) %>%
  ggplot(aes(x = age, y = expected_avg, color = domain)) + 
  geom_smooth(aes(group = interaction(parcellation, domain), size = label), method = "lm", formula = y ~ poly(x, 2), se = FALSE, show.legend = FALSE) +
  geom_smooth(method = "lm", formula = y ~ poly(x, 2), size = 1.25, se = FALSE) +
  scale_color_manual(name = "", values = pal_social_academic) +
  scale_size_manual(values = c(.05, .1, .1)) +
  scale_x_continuous(breaks = c(10, 13, 16)) + 
  scale_y_continuous(breaks = c(-1, 0, 1)) + 
  coord_cartesian(ylim = c(-1.2, 1.2)) +
  facet_grid(~label, labeller = parcel_labeller) +
  labs(x = "\nage", y = "mean predicted BOLD signal value\n") + 
  dcbw +
  theme(legend.position = c(.85, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(0, "cm"))

target_parc_plot = neuro_plot_data %>%
  distinct(parcellation, target, domain, age, label, .keep_all = T) %>%
  group_by(subjectID, age, label, target, parcellation) %>%
  mutate(expected_avg = mean(expected, na.rm = TRUE)) %>%
  ggplot(aes(x = age, y = expected_avg, color = target)) + 
  geom_smooth(aes(group = interaction(parcellation, target), size = label), method = "lm", formula = y ~ poly(x, 2), se = FALSE, show.legend = FALSE) +
  geom_smooth(method = "lm", formula = y ~ poly(x, 2), size = 1.25, se = FALSE) +
  scale_color_manual(values = pal_self_other) +
  scale_size_manual(values = c(.05, .1, .1)) +
  scale_x_continuous(breaks = c(10, 13, 16)) + 
  scale_y_continuous(breaks = c(-1, 0, 1)) + 
  coord_cartesian(ylim = c(-1.2, 1.2)) +
  facet_grid(~label, labeller = parcel_labeller) +
  labs(x = "\nage", y = "mean predicted BOLD signal value\n", color = "") + 
  dcbw +
  theme(legend.position = c(.85, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(0, "cm"))

(h0_fitted = cowplot::plot_grid(domain_parc_plot, target_parc_plot,
                                labels = c('A', 'B'), ncol = 2,
                                rel_widths = c(1, 1)))

hypothesis 1

domain_plot = neuro_plot_data %>%
  group_by(subjectID, age, label, domain, parcellation) %>%
  mutate(expected_avg = mean(expected_mean, na.rm = TRUE)) %>%
  distinct(parcellation, target, domain, age, label, .keep_all = T) %>%
  ggplot(aes(x = age, y = expected_avg, color = domain)) +
  geom_rect(data = subset(label_df, label == "social"), aes(fill = label), color = NA, alpha = .07,
            xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = Inf, show.legend = FALSE) +
  geom_smooth(method = 'lm', formula = y ~ poly(x,2), alpha = .2, se = FALSE, size = 1.25) + 
  scale_y_continuous(breaks = seq(-.2, .45, .2)) +
  coord_cartesian(ylim = c(-.25, .45)) +
  scale_color_manual(values = pal_social_academic) +
  scale_fill_manual(values = "lightgrey") +
  scale_x_continuous(breaks = c(10, 13, 16)) + 
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\nage", y = "mean predicted BOLD signal value\n",  color = '') +
  dcbw +
  theme(legend.position = c(.85, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(0, "cm"))

soc_acad_plot = neuro_plot_data %>%
  group_by(subjectID, age, label, domain, parcellation) %>%
  mutate(expected_avg = mean(expected, na.rm = TRUE)) %>%
  select(subjectID, age, label, domain, expected_avg) %>%
  unique() %>%
  spread(domain, expected_avg) %>%
  mutate(expected_diff = social - academic) %>%
  distinct(parcellation, age, label, .keep_all = T) %>%
  ggplot(aes(x = age, y = expected_diff)) +
  geom_rect(data = subset(label_df, label == "social"), aes(fill = label), alpha = .07,
            xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = Inf, show.legend = FALSE) +
  geom_smooth(aes(group = parcellation, size = label), method = "lm", formula = y ~ poly(x, 2),
              se = FALSE, color = "grey50") +
  geom_smooth(method = 'lm', formula = y ~ poly(x,2), 
              se = FALSE, color = pal_social_academic[2], size = 1.5) + 
  scale_fill_manual(values = "lightgrey") +
  scale_size_manual(values = c(.03, .1, .1)) +
  scale_y_continuous(breaks = seq(-.2, .45, .2)) +
  coord_cartesian(ylim = c(-.25, .45)) +
  scale_x_continuous(breaks = c(10, 13, 16)) + 
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\nage", y = "mean predicted BOLD signal value\n",  color = '') +
  dcbw +
  theme(legend.position = "none")

(h1_fitted = cowplot::plot_grid(domain_plot, soc_acad_plot,
                                labels = c('A', 'B'), ncol = 2,
                                rel_widths = c(1, 1)))

hypothesis 2

target_plot = neuro_plot_data %>%
  group_by(subjectID, age, label, target, parcellation) %>%
  mutate(expected_avg = mean(expected_mean, na.rm = TRUE)) %>%
  distinct(parcellation, target, target, age, label, .keep_all = T) %>%
  ggplot(aes(x = age, y = expected_avg, color = target)) +
  geom_rect(data = subset(label_df, label == "self"), aes(fill = label), color = NA, alpha = .07,
            xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = Inf, show.legend = FALSE) +
  geom_smooth(method = 'lm', formula = y ~ poly(x,2), alpha = .2, se = FALSE, size = 1.25) + 
  scale_y_continuous(breaks = seq(-.2, .3, .1)) +
  coord_cartesian(ylim = c(-.2, .35)) +
  scale_color_manual(values = pal_self_other) +
  scale_fill_manual(values = "lightgrey") +
  scale_x_continuous(breaks = c(10, 13, 16)) + 
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\nage", y = "mean predicted BOLD signal value\n",  color = '') +
  dcbw +
  theme(legend.position = c(.85, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(0, "cm"))

self_other_plot = neuro_plot_data %>%
  group_by(subjectID, age, label, target, parcellation) %>%
  mutate(expected_avg = mean(expected, na.rm = TRUE)) %>%
  select(subjectID, age, label, target, expected_avg) %>%
  unique() %>%
  spread(target, expected_avg) %>%
  mutate(expected_diff = self - other) %>%
  distinct(parcellation, age, label, .keep_all = T) %>%
  ggplot(aes(x = age, y = expected_diff)) +
  geom_rect(data = subset(label_df, label == "self"), aes(fill = label), alpha = .07,
            xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = Inf, show.legend = FALSE) +
  geom_smooth(aes(group = parcellation, size = label), method = "lm", formula = y ~ poly(x, 2),
              se = FALSE, color = "grey50") +
  geom_smooth(method = 'lm', formula = y ~ poly(x,2), 
              se = FALSE, color = pal_self_other[2], size = 1.5) + 
  scale_fill_manual(values = "lightgrey") +
  scale_size_manual(values = c(.03, .1, .1)) +
  scale_y_continuous(breaks = seq(-.2, .3, .1)) +
  coord_cartesian(ylim = c(-.2, .35)) +
  scale_x_continuous(breaks = c(10, 13, 16)) + 
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\nage", y = "mean predicted BOLD signal value\n",  color = '') +
  dcbw +
  theme(legend.position = "none")

(h2_fitted = cowplot::plot_grid(target_plot, self_other_plot,
                                labels = c('A', 'B'), ncol = 2,
                                rel_widths = c(1, 1)))

hypothesis 3

int_plot = neuro_plot_data %>%
  distinct(parcellation, target, domain, age, label, .keep_all = T) %>%
  ggplot(aes(x = age, y = expected_mean, group = interaction(target, domain), color = domain, linetype = target)) +
  geom_smooth(method = 'lm', formula = y ~ poly(x,2), alpha = 1, se = FALSE, size = 1.25) + 
  scale_y_continuous(breaks = c(-.4, -.2, 0, .2, .4)) +
  coord_cartesian(ylim = c(-.4, .5)) +
  scale_color_manual(values = pal_social_academic) +
  scale_linetype_manual(name = "", values = c("dotted", "solid")) +
  scale_x_continuous(breaks = c(10, 13, 16)) + 
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\nage", y = "mean predicted BOLD signal value\n",  color = '', linetype = '') +
  guides(linetype = guide_legend(override.aes = list(color = "black"))) +
  dcbw +
  theme(legend.position = c(.75, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(.5, "cm"),
        legend.direction = "vertical",
        legend.box = "horizontal")

soc_acad_int_plot = neuro_plot_data %>%
  group_by(subjectID, age, label, target, domain, parcellation) %>%
  mutate(expected_avg = mean(expected, na.rm = TRUE)) %>%
  select(subjectID, age, label, target, domain, expected_avg) %>%
  unique() %>%
  spread(domain, expected_avg) %>%
  mutate(expected_diff = social - academic) %>%
  distinct(parcellation, age, label, .keep_all = T) %>%
  ggplot(aes(x = age, y = expected_diff, color = target)) +
  geom_smooth(aes(group = interaction(parcellation, target), size = label),
              method = "lm", formula = y ~ poly(x, 2), se = FALSE, show.legend = FALSE) +
  geom_smooth(method = 'lm', formula = y ~ poly(x,2), 
              se = FALSE, size = 1.5) + 
  scale_color_manual(values = pal_self_other) +
  scale_size_manual(values = c(.03, .1, .1)) +
  scale_y_continuous(breaks = c(-.4, -.2, 0, .2, .4)) +
  coord_cartesian(ylim = c(-.4, .5)) +
  scale_x_continuous(breaks = c(10, 13, 16)) + 
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\nage", y = "mean predicted BOLD signal value\n",  color = '') +
  dcbw +
  theme(legend.position = c(.85, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(0, "cm"))


self_other_int_plot = neuro_plot_data %>%
  group_by(subjectID, age, label, domain, target, parcellation) %>%
  mutate(expected_avg = mean(expected, na.rm = TRUE)) %>%
  select(subjectID, age, label, domain, target, expected_avg) %>%
  unique() %>%
  spread(target, expected_avg) %>%
  mutate(expected_diff = self - other) %>%
  distinct(parcellation, age, label, .keep_all = T) %>%
  ggplot(aes(x = age, y = expected_diff, color = domain)) +
  geom_smooth(aes(group = interaction(parcellation, domain), size = label),
              method = "lm", formula = y ~ poly(x, 2), se = FALSE, show.legend = FALSE) +
  geom_smooth(method = 'lm', formula = y ~ poly(x,2), 
              se = FALSE, size = 1.5) + 
  scale_color_manual(values = pal_social_academic) +
  scale_size_manual(values = c(.03, .1, .1)) +
  scale_y_continuous(breaks = c(-.4, -.2, 0, .2, .4)) +
  coord_cartesian(ylim = c(-.4, .5)) +
  scale_x_continuous(breaks = c(10, 13, 16)) + 
  facet_grid(~ label, labeller = parcel_labeller) +
  labs(x = "\nage", y = "mean predicted BOLD signal value\n",  color = '') +
  dcbw +
  theme(legend.position = c(.85, .15),
        legend.spacing.y = unit(.01, 'cm'),
        legend.margin = unit(0, "cm"))

(h3_fitted = cowplot::plot_grid(int_plot, soc_acad_int_plot, self_other_int_plot,
                       labels = c('A', 'B', 'C'), ncol = 3))

compare raw and fitted models

hypothesis 0

raw

h0_raw

fitted

h0_fitted

hypothesis 1

raw

h1_raw

fitted

h1_fitted

hypothesis 2

raw

h2_raw

fitted

h2_fitted

hypothesis 3

raw

h3_raw

fitted

h3_fitted